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ABSTRACT 

We model the distribution of neutral Hydrogen (HI hereafter) in the post-reionization universe. 
This model uses gravity only N-Body simulations and an ansatz to assign HI to dark matter 
haloes that is consistent with observational constraints and theoretical models. We resolve 
the smallest haloes that are likely to host HI in the simulations, care is also taken to ensure 
that any errors due to the finite size of the simulation box are small. We then compute the 
smoothed one point probability distribution function and the power spectrum of fluctuations 
in HI. This is compared with other predictions that have been made using different techniques. 
We highlight the significantly high bias for the HI distribution at small scales. This aspect has 
not been discussed before. We then discuss the prospects for detection with the MWA, GMRT 
and the hypothetical MWA5000. The MWA5000 can detect visibility correlations at large 
angular scales at all redshifts in the post-reionization era. The GMRT can detect visibility 
correlations at lower redshifts, specifically there is a strong case for a survey at z ~ 1.3. We 
also discuss prospects for direct detection of rare peaks in the HI distribution using the GMRT. 
We show that direct detection should be possible with an integration time that is comparable 
to, or even less than, the time required for a statistical detection. Specifically, it is possible to 
make a statistical detection of the HI distribution by measuring the visibility correlation, and, 
direct detection of rare peaks in the HI distribution at z ~ 1.3 with the GMRT in less than 
1000 hours of observations. 

Key words: methods: N-Body simulations, cosmology: large scale structure of the universe, 
galaxies: evolution, radio-lines: galaxies 



1 INTRODUCTION 

Large scale structures in the universe are believed to have formed 
by gravitational amplification of small perturbations (Peebles 1980; 
Shandarin & Zeldovich 1989; Peacock 1999; Padmanabhan 2002; 
Bernardeau et al. 2002). Much of the matter in galaxies and clus- 
ters of galaxies is the so called dark matter that is believed to be 
weakly interacting and non-relativistic (Trimble 1987; Komatsu et 
al. 2009). Dark matter responds mainly to gravitational forces, and 
by virtue of larger density than baryonic matter, assembly of mat- 
ter into haloes and large scale structure is driven by gravitational 
instability of initial perturbations. 

Galaxies are believed to form when gas in highly over-dense 
haloes cools and collapses to form stars in significant numbers 
(Hoyle 1953; Rees & Ostriker 1977; Silk 1977; Binney 1977). 
The formation of first stars (McKee & Ostriker 2007; Zinnecker 
& Yorke 2007; Bromm & Larson 2004) in turn leads to emission of 
UV radiation that starts to ionize the inter-galactic medium (IGM). 
The period of transition of the IGM from a completely neutral to 



a completely ionized state is known as the epoch of reionization 
(EoR), e.g., see Loeb & Barkana (2001). Observations indicate that 
the process of reionization was completed before z ~ 6 (Fan, Car- 
illi, & Keating 2006; Becker et al. 2001; Fan et al. 2006). Many pos- 
sible sources of ionizing radiation have been considered, although 
stellar sources are believed to be the most plausible candidates, see, 
e.g., Bagla, Kulkarni, & Padmanabhan (2009). 

Prior to the EoR, almost all the Hydrogen in the universe is 
in atomic form. Through the EoR Hydrogen is ionized till we are 
left with almost no HI in the inter-galactic medium (IGM) and al- 
most all the HI resides in the inter stellar medium (ISM) of galax- 
ies. We focus on the post-reionization era in this work and our aim 
is to make predictions about the distribution of HI in this regime. 
It had been proposed by Sunyaev & Zeldovich (1972, 1975) that 
the hyperfine transition of Hydrogen may be used to probe primor- 
dial galaxies. This problem has been approached in past from the 
perspective of making specific predictions for existing instruments 
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like the Giant Meterwave Radio Telescope (GMRT) 1 (Subrama- 
nian & Padmanabhan 1993; Kumar, Padmanabhan, & Subramanian 
1995; Bagla, Nath, & Padmanabhan 1997; Bagla 1999; Bharadwaj, 
Nath, & Sethi 2001 ; Bharadwaj & Sethi 2001 ; Bagla & White 2003; 
Bharadwaj & Srikant 2004; Bharadwaj & Ali 2005), or upcoming 
instruments like the the MWA 2 , ASKAP a MeerKAT 4 and SKA 5 , 
e.g., see Wyithe, Loeb, & Geil (2008). Much work in the last decade 
has focused on making predictions for the power spectrum of fluc- 
tuations in HI, with an implicit assumption that it is easier to make 
a statistical detection than a direct detection. Further, it has been 
argued that HI can be used as a tracer of the large scale structure 
and observations can be used to constrain cosmological parameters 
with a special emphasis on observations of the baryon acoustic os- 
cillations (BAO) in the power spectrum (Visbal, Loeb, & Wyithe 
2008; Chang et al. 2008; Bharadwaj, Sethi, & Saini 2009). Most of 
the upcoming instruments are sensitive to the power spectrum of 
fluctuations at large scales, certainly larger than the scale of non- 
linearity at the relevant epoch, and hence a significant fraction of 
the work done in terms of making predictions is based on linear 
theory or approximations that work well in the linear and quasi- 
linear regime. The halo model has also been used for predicting the 
spectrum of fluctuations in the post-reionization universe (Wyithe 
et al. 2009; Wyithe & Brown 2009). 

In this study we revisit the issue of predicting fluctuations and 
use high resolution N-Body simulations. This allows us to study 
fluctuations at scales comparable to, or even smaller than the scale 
of non-linearity. We use dark matter simulations along with simple 
ansatz for assigning neutral Hydrogen in order to make predictions 
of fluctuations in surface brightness temperature. Our model is dis- 
cussed in detail in §2 where we summarize our knowledge of the HI 
distribution at high redshifts and motivate the assignment schemes 
we use in this work . The N-Body simulations used by us are de- 
scribed in §2. 1 . Results of the HI signal are presented in §3. We then 
move on to review the relation between the source flux and observ- 
ables in radio interferometers, i.e. visibilities, in §4 and also discuss 
the sensitivity of interferometers. Finally we look at the prospects 
of detection, both statistical and of rare peaks in §5. We conclude 
with a discussion in §6. 



2 MODELING THE HI DISTRIBUTION 

In this section we describe our model of the HI distribution at high 
redshifts. Our knowledge of the HI distribution in the universe is 
derived mainly from QSO absorption spectra, where the gas ab- 
sorbs in the Lyman-a transition of the Hydrogen atom. We know 
from observations of these absorption spectra that much of the 
inter-galactic medium (IGM) is highly ionized and does not contain 
a significant amount of neutral Hydrogen. Most of the neutral Hy- 
drogen resides in relatively rare damped Lyman-a systems (DLAS) 
(Wolfe, Gawiser, & Prochaska 2005). DLAS and other high column 
density absorption features are believed to arise due to gas within 



1 See http://gmrt.ncra.tifr.res.in/ for further details. 

2 The Murchison Widefield Array (MWA). More details are available at 
http://www.mwatelescope.org/ 

3 The Australian Square Kilometer Array Pathfinder. See 
http://www.atnf.csiro.au/projects/askap/ for details. 

4 The South African Square Kilometer Array Pathfinder. See 
http://www.ska.ac.za/meerkat/overview.php for details. 

5 The Square Kilometer Array (SKA). See http://www.skatelescope.org/ 
for details. 



galaxies (Haehnelt, Steinmetz, & Rauch 2000; Gardner et al. 2001). 
It is possible to make a quantitative estimate of the total neutral Hy- 
drogen content in DLAS and study the evolution of the total neutral 
Hydrogen content of the universe (Storrie-Lombardi, McMahon, & 
Irwin 1996; Rao & Turnshek 2000; Peroux et al. 2005). These ob- 
servations indicate that at 1 < z < 5, the neutral Hydrogen con- 
tent of the universe is almost constant with a density parameter of 
Qhi ^ 0.001 6 . 

At low redshifts, the HI content can be estimated more di- 
rectly through emission in the Hyperfine transition. Observations 
in the local universe indicate a much lower neutral Hydrogen con- 
tent than seen at z > 1 (Zwaan et al. 2005). The neutral gas frac- 
tion in galaxies at intermediate redshifts appears to be much higher 
than in galaxies in the local universe (Lah et al. 2007, 2009), i.e., 
the neutral gas fraction in galaxies appears to increase rapidly with 
redshift z. 

The spin temperature couples to the gas temperature through 
collisions of atoms with other atoms, electrons, ions and also the 
Wouthuysen-Field effect (Purcell & Field 1956; Field 1958, 1959; 
Wouthuysen 1952a,b; Zygelman 2005; Furlanetto, Oh, & Briggs 
2006; Furlanetto & Furlanetto 2007a,b). Observations of 21 cm ab- 
sorption by DLAS indicate that the spin temperature is orders of 
magnitude higher than the temperature of the cosmic microwave 
background radiation (CMBR) at corresponding redshifts (Chen- 
galur & Kanekar 2000; Kanekar et al. 2009). This implies that the 
emission in the 21 cm hyperfine transition can be safely assumed 
to be proportional to the density of neutral Hydrogen, see, e.g., 
Furlanetto, Oh, & Briggs (2006). 
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where xhi is the fraction of Hydrogen in neutral form, 8 is the 
density contrast of the gas distribution, T cm b is the temperature of 
the CMBR T s is the spin temperature defined using the relative 
occupation of the two levels for the hyperfine transition: 
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where subscripts 1 and correspond to the excited and ground state 
levels of the hyperfine transition, = hv/ks = 68 mK is the 
temperature corresponding to the transition energy, and (<?i/<?o) = 
3 is the ratio of the spin degeneracy factors of the levels. 

Observations indicate that neutral gas is found only in galaxies 
in the post reionization universe. We also know that at very low red- 
shifts galaxies in groups and clusters do not contain much neutral 
gas. As cold gas is associated with galaxies, we may assume that 
it exists only in haloes that are more massive than the Jeans mass. 
Further, we may assume that neutral gas is predominantly found 
in galaxies and not in larger haloes that may contain several large 
galaxies. Jeans mass for haloes in a photo-ionizing UV background 



6 These observational constraints are given as the ratio of the comoving 
density of neutral Hydrogen to the present day critical density. We adopt 
the same convention here. 
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depends on the shape of the spectrum of the ionizing radiation, typ- 
ically we expect gas in haloes with a circular velocity in excess of 
60 km/s to cool, fragment and form stars. The mass within the virial 
radius is related to the circular velocity and the collapse redshift as: 

Simulation studies show that DLAS can reside in haloes with even 
lower circular velocities (Pontzen et al. 2008), i.e., haloes with mass 
lower than the expected Jeans mass can contain significant amount 
of neutral Hydrogen. The gas in these haloes is able to self shield 
and maintain a significant amount of HI even though the amount 
of gas is insufficient for sustaining star formation 7 . We use this in- 
put and impose a lower cutoff in circular velocity of 30 km/s, i.e., 
haloes with a lower circular velocity are not assigned any HI. This 
lower limit is appropriate for z ~ 3 and it is likely that the limit 
shifts to higher values at lower redshifts. We keep the value fixed 
at 30 km/s in order to make a conservative estimate of clustering of 
such haloes. Pontzen et al. (2008) also found that haloes much more 
massive than a few times 10 11 M Q do not host significant amount 
of HI. This is consistent with observations in the local universe 
where galaxies in groups and clusters of galaxies contain very little 
HI. Some observations suggest that the neutral fraction in galaxies 
in outer parts of clusters of galaxies at intermediate redshifts may 
be high compared to corresponding galaxies in the local universe 
(Lah et al. 2009). Assignment for HI in more massive haloes is 
done in such a way as to ensure that very massive haloes have a 
zero or negligible fraction of the total gas in neutral form. In case 
larger haloes have more HI than we assume, our results under esti- 
mate the true signal and may be considered as a lower bound. The 
transition scale for the higher masses is chosen to coincide with 
200 km/s. Given that it is easier to estimate mass than circular ve- 
locity for haloes identified in simulations, in particular for haloes 
with a small number of particles, we choose to translate the thresh- 
old in circular velocity to a threshold in mass assuming the collapse 
redshift and the redshift at which the system is being observed to be 
the same, i.e., z c = z. While this may lead to inaccuracies regard- 
ing inclusion of haloes close to the low mass end, we do not expect 
it to influence the results in a significant manner. The minimum and 
the maximum masses for the redshifts considered here are listed in 
Table 1. 

Given a halo of mass M, all particles in it are assigned an 
equal HI mass which is a fraction of their total mass. We describe 
three kinds of mass assignments here. 

Fl(M) = /l (M m i n < M < Mm ax ) (4) 

F 2 (M) = -fi — 2 (M min <M) (5) 

1+ ( M \ 

F 3 (M) = —fr ^ (M mm <M) (6) 



7 Haloes that collapse much earlier may have hosted star formation, in- 
deed the ISM of such galaxies may have been blown away due to feedback. 
However, few haloes with a low circular velocity that collapse very early are 
likely to survive without merger well into the post-reionization era. Errors 
arising out of HI assignment to low circular velocity haloes should not lead 
to significant changes in our calculations. Errors arising from this will gen- 
erally lead to an under estimate of the signal (see discussion of dependence 
of bias on the lower mass limit), hence our predictions may be considered 
to be a lower limit on the signal. 



Table 1. The table lists the minimum mass of haloes M min that hosts HI 
in the simulated HI distribution. The first column lists the redshift where 
simulated maps are constructed, the second column lists M min we get by 
assuming a lower bound of 30 km/s on the circular velocity of haloes, and 
column 4 lists the actual minimum halo mass used in construction of the 
simulated maps. All masses are given in units of solar mass Mq. Column 
3 lists the characteristic maximum mass of haloes computed using an upper 
bound on circular velocity of 200 km/s and Column 5 shows the actual 
upper bound used while making simulated maps. Column 6 lists the HI 
fraction of the total baryonic mass if we use the first HI assignment scheme 
(see Eqn. 4. 
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Here F k (M) is the mass fraction of HI in a given halo. The frac- 
tion of baryons that is in the form of neutral Hydrogen is given by 
Fk{M)Q, nr /Q.b, where fib is the density parameter of Baryons and 
Q„ r is the density parameter for non-relativistic matter. For the best 
fit model for WMAP-5 (Komatsu et al. 2009) the ratio has the value 
Q nr /D. b ~ 1/0.17 ~ 6. 

In all three cases the constant ft is determined by normaliz- 
ing the HI mass in the simulation volume to SIhi = 0.001, the 
density parameter of HI indicated by observations at relevant red- 
shifts (Peroux et al. 2005). It is noteworthy that this normalization 
requires a significant fraction of gas in haloes in the mass range 
between M m ; n and M max be neutral at high redshifts. In the as- 
signment scheme with the sharp cutoff (scheme 1, Eqn. 4), the HI 
fraction is the same in all haloes with mass between M min and 
M max . In the second assignment scheme the HI fraction in a halo 
decreases monotonically with increasing halo mass, and for large 
halo masses the HI mass goes to zero. This scheme applies to the 
physical situation where very large haloes do not have any neutral 
Hydrogen. In scheme 3, the HI fraction decreases monotonically 
with increasing halo mass but the HI mass goes to a constant value. 
The motivation of this scheme is to allow for some neutral Hy- 
drogen surviving in massive haloes in galaxies. In the following 
section, we compare the resulting distribution of HI with the three 
assignment schemes described above. The present work is the first 
one where mass resolution of simulations is adequate for resolving 
the smallest haloes that may host HI and the simulation volume is 
also sufficient to make statements about the large scale distribution. 
We also ensure that the finite size of the simulation volume does not 
affect the results presented here. 



2.1 AT -Body Simulations 

We use gravity only simulations run with the TreePM code (Bagla 
2002; Bagla & Ray 2003; Khandai & Bagla 2009). The suite of 
simulations used here is described in the Table 2. The cosmological 
model and the power spectrum of fluctuations corresponds to the 
best fit model for WMAP-5: Q. nr = 0.26, Cl A = 0.74, n s = 0.96, 
(T 8 = 0.79, h = 0.72, n b h 2 = 0.02273 (Komatsu et al. 2009). 

We use the Friends-of-Friends (FOF) (Davis et al. 1985) algo- 
rithm with a linking length I = 0.2 to identify haloes and construct 
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Table 2. Columns 1 and 2 list the size of the box and the number of particles 
used in the simulations. Columns 3 and 4 give the mass and force resolution 
of the simulations, while columns 5 and 6 tell us the redshift at which the 
simulations were terminated and the redshifts for which the analysis were 
done 
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a halo catalog. The HI assignment schemes are then used to obtain 
the distribution of neutral Hydrogen. 

Several existing and upcoming instruments can probe the post- 
reionization universe using redshifted HI emission. The GMRT can 
observe redshifted HI emission from a few selected redshift win- 
dows whereas most other instruments have continuous coverage in 
redshift bounded on two sides. We choose to focus on the GMRT 
windows, as these are representative of the range of redshifts in 
the post-reionization universe. In particular we will focus on the 
following redshift windows of GMRT: z out = 5.04, 3.34, 1.33. 
Finite box effects can lead to significant errors in the distribution 
of haloes that host galaxies, apart from errors in the abundance of 
haloes of different masses (Bagla & Prasad 2006; Bagla, Prasad, 
& Khandai 2009). The choice of simulations used in this work 
ensures that such effects do not contribute significantly. Previous 
studies have indicated that at z ~ 0, we need a simulation box 
with Lbox > 140 h _1 Mpc for the finite size effects to be negligi- 
ble (Bagla & Ray 2005; Bagla & Prasad 2006). On the other hand 
the mass resolution of particles decreases as the cube of simulation 
volume. We balance the requirements of high mass resolution and a 
sufficiently large box size by using different simulations for study- 
ing the HI distribution at different redshifts. Details of the simula- 
tions are given in the table 2. 



3 RESULTS 

Given an HI mass assigment for particles in a simulation, we can 
proceed to compute the expected signal from the HI distribution by 
making mock radio maps and spectra. We also compute the power 
spectrum in both the real and the redshift space. For redshift space 
calculations, we use the peculiar velocity information of particles 
in haloes. Unlike earlier studies, we resolve haloes of individual 
galaxies and hence the internal velocity dispersion is naturally ac- 
counted for and there is no need to add it by hand (Kumar, Padman- 
abhan, & Subramanian 1995; Bagla, Nath, & Padmanabhan 1997; 
Bagla & White 2003). 

We use the clouds-in-cell (CIC) smoothing to interpolate den- 
sities from particle locations to grid points for the purpose of com- 
puting densities and the power spectrum. It is convenient to ex- 
press the HI power spectrum in terms of the brightness tempera- 
ture, though in comparisons with the dark matter power spectrum 
we revert to the usual dimensionless form. 

We define the real and redshift space scale dependant bias by 
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Figure 1. Effect of HI mass assignment type on bias. Solid, dashed and dot- 
dashed lines are for HI assignment types 1, 2 and 3 as described in eqs. 4-6. 
Bias is computed at z = 3.34 for the 51.2h — *Mpc box with M max = 
W 115 M . . 
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We start by checking the effect of HI assignment scheme on the 
resulting distribution. We have computed the bias b(k) at z = 3.34 
for the Li,ox = 51.2h _1 Mpc simulation using the three schemes 
introduced above. The results are shown in Figure 1 for the three 
assignment types described in Eqs. (4-6). We see that at large scales 
the three assignment schemes give very similar results, whereas 
there is some disagreement between the assigment scheme Eqn.(4) 
and the other two. The difference can be attributed to the fact that 
the scheme described in Eqn.(4) puts more HI mass in haloes with 
masses near M max . The differences between the HI assignment 
schemes are relatively minor. We choose to work with the scheme 
described in Eqn.(6) due to a better physical justification, as dis- 
cussed in §2. As mentioned before, if large mass haloes happen to 
have a larger HI fraction than we assume then the signal will be 
larger than what we predict here. 

It is noteworthy that while the bias is strongly scale dependent 
at large k (small scales), it flattens out to a constant value at small 
k (large scales). 

In order to check for the effects of a finite box size, we car- 
ried out a test for the 23 h _1 Mpc simulation box. Instead of using 
the range of halo masses for HI assignment that is appropriate for 
z — 5.1, we work with a slightly smaller range so that haloes of 
these masses can be found in the 51.2 h _1 Mpc simulation as well. 
Figure 2 shows the dark matter power spectrum (Top panel) and the 
HI power spectrum (Lower panel) with the HI assignment restricted 
to a smaller range of masses, as described above. We see that the 
dark matter power spectra from the two simulations agree through 
the range of scales where there is an overlap. The HI power spec- 
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Figure 2. Top: The dark matter power spectrum in the simulations with 
Lbox = 23 h _1 Mpc and L^ ox = 51.2 h — 1 Mpc at z = 5.1. The agree- 
ment in the two curves indicates that the finite box size effects do not lead to 
errors in the smaller box at the redshift of interest. Bottom: HI power spec- 
trum in the same simulations, see text for details. Show the effect of finite 
box size on the simulated HI power spectrum, good agreement between the 
two curves imples that the error is small. 



tra also agree, though not as well as the dark matter power spec- 
tra. These differences are so small that we do not expect these to 
affect the final results in a significant manner. The difference can 
be attributed to the fact that clustering of haloes is affected more 
strongly by the box size effects in simulations (Gelb & Bertschinger 
1994; Bagla & Ray 2005). 

These tests validate our approach for assignment of HI to 
haloes, and also show that the effects of a finite box-size are not 
significant at the level of the power spectrum or the mass function. 



o _ 
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Figure 3. The power spectrum of fluctuations: the solid line shows the lin- 
early extrapolated power spectrum, the dashed line shows the non-linear 
dark matter power spectrum and the dot-dashed line shows the HI power 
spectrum. All power spectra are for z = 3.34. The dark matter and the 
HI power spectra have been computed with the simulation with Lt, ox = 
Sl^h^Mpc. 




k(hMpc -1 ) 

Figure 4. The dependance of bias on Af m ; n and M max . The solid curve 
shows b(k) for the values of M min and M max shown in Table 1. The 
dashed line shows b(k) when we use the reference value for M max but 
increase M m i n to twice the reference value. The dot-dashed curve shows 
b(k) when M max is chosen to be higher: 10 Mq, while M m j n is kept 
at the reference alue. We see that as the characteristic mass of haloes with 
HI increases, b(k) increases. All curves are for z = 3.34 and have been 
computed with the simulation with Lt, ox = 51.2 h _1 Mpc. 
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Figure 5. This figure shows a scatter plot of Shi smoothed at a scale of 
3h~ ijvlpc with spherical top hat window, plotted as a function of &dm 
smoothed at the same scale. The figure shows a random subset of points 
from the simulation with L^ox 

= Sl^h^Mpcatz = 3.34. 

3.1 Bias 

This is amongst the first studies of the HI distribution at high red- 
shifts where we resolve the smallest haloes that can host significant 
amount of HI while ensuring that the finite box size effects do not 
lead to an erronous distribution of haloes. One of the points that 
we can address here is the effect of non-linear clustering on the HI 
distribution and the scale dependance of bias. Some of these effects 
are illustrated in Figure 3 where we have plotted the power spec- 
trum of fluctuations: the solid line shows the linearly extrapolated 
power spectrum, the dashed line shows the non-linear dark matter 
power spectrum and the dot-dashed line shows the HI power spec- 
trum. All power spectra are for z = 3.34. The dark matter and 
the HI power spectra have been computed with the simulation with 
Lbox = 51.2 h _1 Mpc. 

It is apparent that at k > 0.5 h Mpc -1 the effects of non-linear 
clustering significantly enhance the dark matter power spectrum. At 
k ~ 10 h Mpc^ 1 , the enhancement is close to an order of magni- 
tude. This is of utmost interest for upcoming instruments that can 
resolve small angular scales. 

We find that bias b(k) for the HI distribution is much greater 
than unity at high redshifts. This is to be expected of galaxies at 
high redshifts (Fry 1996; Mo & White 1996; Bagla 1998a,b; Mo, 
Mao, & White 1999; Baugh et al. 1999; Magliocchetti et al. 2000; 
Benson et al. 2000; Roukema & Valls-Gabaud 2000; Sheth, Mo, 
& Tormen 2001; Wyithe & Brown 2009). We note that the bias is 
scale dependent, and leads to a larger enhancement in the power 
spectrum at very small scales. 

The value of bias depends strongly on the choice of the char- 
acteristic mass of haloes with HI. This is shown in Figure 4. All 
curves are for z = 3.34 and have been computed with the simula- 
tion with Lbox = 51.2 h _1 Mpc. This figure illustrates that the bias 
at all scales varies monotonically with the charactersitic mass of 
haloes with HI. Variation is gentle at large scales but fairly strong 
at small scales. Therefore it is extremely important to have an accu- 



rate estimate of the charactersitic masses of such haloes. Indeed, it 
has been pointed out that observations of the amplitude of cluster- 
ing in the HI distribution can be used to constrain masses of haloes 
that host DLAS (Wyithe 2008). 

While the preceeding figures describe the statistical bias com- 
puted from the ratio of power spectra, Figure 5 shows the stochas- 
ticity of bias (Dekel & Lahav 1999) in the HI distribution. This 
figure shows a scatter plot of Shi smoothed at a scale of 3 h _1 Mpc 
with a spherical top hat window, plotted as a function of Sdm 
smoothed at the same scale. The figure shows a random subset of 
points from the simulation with Lbox = 51.2 h _1 Mpc at z = 3.34. 
The scatter about the average trend in the Shi — Sdm is significant, 
and increases as we go towards large overdensities in dark matter. 
The scatter becomes small if we smooth the density distribution at 
much larger scales. 

Similar results concerning the effect of non-linear gravita- 
tional clustering and bias are obtained for other redshifts. We con- 
centrate on the evolution of bias and the power spectrum rather than 
go into a very detailed discussion of the HI distribution at each red- 
shift. We discuss simulated maps in the following subsection. 

Figure 6 (left panel) shows the evolution of the redshift space 
power spectrum for the HI distribution. Curves show the A H i(k) 
as a function of k for z = 5.1 (solid line), z = 3.34 (dashed 
line) and z — 1.3 (dot-dashed line). Unlike the power spectrum in 
real space, the enhancement at very small scales is less strong and 
this is due to velocity dispersion within haloes at smaller scales. 
Power at larger scales is enhanced due to the Kaiser effect (Kaiser 
1987), leading to a relatively stronger enhancement at larger scales, 
hence the difference in shapes of the HI power spectrum in the 
real space and the redshift space. Even though the amplitude of 
density perturbations increases with time, we see that the amplitude 
of brightness temperature fluctuations decreases instead. This is a 
manifestation of the decreasing correlation function for galaxies, 
see, e.g., Bagla (1998b); Roukema & Valls-Gabaud (1998). 

The evolution of bias at large scales, denoted as bu n (k) to em- 
phasise that it refers to scales where clustering is still in the linear 
regime, is shown in the right panel of Figure 6. We see that this de- 
creases from around 2.5 at z = 5.1 to 1.2 at z = 1.3, the variation 
is slightly steeper than (1 + z) 1/2 and hence the gradual lowering 
of the amplitude of the brightness temperature power spectrum. 

3.2 Radio Maps 

Figure 7 shows the simulated radio maps and spectra. The usual 
conversion from HI density to signal has been used for this (Furlan- 
etto, Oh, & Briggs 2006). The top row shows simulated radio maps, 
the top-left panel is for z = 3.34 and the top-right panel is for 
z = 1.3. The pixel size is chosen to be the same as the resolution 
of the central square of the GMRT. The bandwidth of the map is 
0.5MHz and 1.0MHz for z = 3.34, 1.33 respectively. The corre- 
sponding contours mark regions with signal (9.2, 5.75, 1.15) ^iJy 
for z = 3.34 and (16.8, 10.5, 2.1) ^Jy for z = 1.33. As we can 
see in these maps, bright regions have a large angular extent. This 
fact may be used to enhance prospects of detection §5. 

The lower panel shows spectra from the mock radio maps. 
These spectra have been plotted for regions that are the size of the 
resolution of the central square of the GMRT. We have chosen to 
plot the spectrum through the brightest region in the simulation. 
The lower-left panel corresponds to z = 3.34 and the lower right 
panel shows the spectrum for z = 1.3. The peak flux is around 
0.025 mjy and the FWHM of the line is close to 300 kHz. The 
strength of the signal for z = 3.34 is comparable to what we found 
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Figure 6. Left: Evolution of the HI powerspectrum (k). Solid, dashed, dot-dashed lines are for redshifts z = 5.04, 3.34, 1.33 respectively. Right: 
Evolution of the HI linear bias. 



in an earlier study, although the line width is much smaller. The 
reason for a smaller line width is likely to be the supression of the 
HI fraction in very massive haloes. Considering the peak for which 
the spectrum has been plotted here, it may be possible to make a 
2 a detection with about 2 x 10 3 hours of observations with the 
central square of the GMRT (§5). Given that the GMRT observes a 
much larger volume at these redshifts that the simulation volume, 
it is highly likely that even rarer peaks in the HI distribution will 
be observed in a generic pointing (Subramanian & Padmanabhan 
1993). 

The strentgh of the predicted signal at z = 1.3 is of the order 
of 0.05 mjy, with an FWHM of around 1 MHz. 

Detection of rare peaks in the HI distribution is an exciting 
possibility. The size of the region represented in these rare peaks 
is fairly large, and it should be possible to establish the total mass 
contained in these regions using observations in other wavebands. 
This will allow us to estimate Q.hi in emission, and hence provide 
an independent measurement of the cold gas fraction. 

The time required for detection of HI at high redshifts with the 
GMRT increases rapidly as we go to higher redshifts. This trend 
continues as we move to z ~ 5.1 and therefore we do not discuss 
those results in detail here. 

Instruments like the GMRT, MWA can in principle detect sig- 
nal from the HI distribution at high redshifts. The angular resolu- 
tion of the MWA is fairly poor and hence the effects of non-linear 
clustering and scale dependant bias do not make a significant im- 
pact on predictions. We shall discuss the prospects for detection 
with the GMRT and the MWA in detail in the §5. 



4 SIGNAL AND NOISE IN INTERFEROMETERS 

In this section we briefly review the relation between the flux from 
sources and the observable quantities for radio interferometers. We 
then proceed to a discussion of the sensitivity of interferometers 
and the corresponding limitations arising from that. 



4.1 Visibility Correlation 

The quantity measured in radio-interferometric observations is the 
visibility V(U, v) which is measured in a number of frequency 
channels v — (y + Av ) across a frequency bandwidth B for every 
pair of antennas in the array. The visibility is related to the sky 
specific intensity pattern I v (&) as 

V(U,u) = J(f&A(&)I u (&)e 2m&u . (9) 

The baseline U is d/A where d is the antenna separation projected 
in the plane perpendicular to the line of sight. is a two dimen- 
sional vector in the plane of the sky with origin at the center of 
the field of view, and A(&) is the beam pattern of the individual 
antenna. For the GMRT this can be well approximated by Gaus- 
sian A(&) = e- e2/t>0 \ where 6 = |0|. We use 6 = 0.56° at 
610 MHz for the GMRT, and it scales as the inverse of frequency. 

We next consider the visibility- visibility correlation (hereafter 
only visibility correlation) measured at two different baselines U 
and U + AU and at two frequencies v and v + Av. As ar- 
gued in Bharadwaj & Ali (2005) the visibilities at baselines U and 
U + AU will be correlated only if AU < l/n0 o . Visibilities 
at different frequencies are expected to be correlated only if the 
flux from sources is correlated. This is certainly true of continuum 
radiation. In case of spectral lines, like the redshifted HI 21 cm 
line being considered here, visibilities are expected to be correlated 
over a range of frequencies comparable to width of spectral lines. 
For generic applications, we can take this to be about 1 MHz though 
we can use the simulated radio maps for a more refined model. As a 
first approximation, the visibilities at frequencies v and v + Av can 
be assumed to be uncorrelated for Av > 1 MHz (Datta, Choud- 
hury, & Bharadwaj 2007) for range of baselines of our interest. 
We use the power spectra to compute visibility correlations using 
the relation enunciated by Bharadwaj & Sethi (2001). To calculate 
the strength and nature of visibility correlation of the HI fluctua- 
tions we consider the visibility correlation at same baselines and 
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Figure 7. The top row shows mock radio maps, the top-left panel is for z = 3.34 and the top-right panel is for z = 1.3. The pixel size is chosen to be the 
resolution of the central square of the GMRT. The bandwidth of the map is 0.5MHz and 1.0MHz for z = 3.34, 1.33 respectively. The corresponding contours 
mark regions with signal (9.2, 5.75 ,1.15) fi]y and (16.8, 10.5, 2.1) /ijy for z = 3.34 and 1.33 respectively. The mock maps show the brightest regions from 
the simulation. As we can see, these regions have a large angular extent this fact may be used to enhance prospects of detection. The lower panel shows spectra 
from the mock radio maps. These spectra have been plotted for regions that are the size of the resolution of the central square of the GMRT for z = 3.34 and 
twice the resolution of the central square of GMRT for z = 1.33. We chose to plot the spectrum through the brightest region in the simulation. The lower-left 
panel corresponds to z = 3.34 and the lower right panel shows the spectrum for z = 1.3. The strength of the signal for z = 3.34 is comparable to what 
we found in an earlier study. The strentgh of the predicted signal at z = 1.3 is or the order of 0.05 mJy, with an FWHM of around 1 MHz. The best rms 
sensitivity achieved in this band at the GMRT is 0.02 mJy. This has been achieved with a bandwidth ~ 32 MHz. 



frequencies which can be written as (for details see Bharadwaj & 
Ali (2005)) 

r2/i2 roo 

vv*(\j, v) = (V(u, v)v*(v, u)) = -0rj o dk n P m (k) 

(10) 

where I v is mean HI specific intensity, r v is comoving distance cor- 
responds to frequency v = 1420/ (1 + z). k = feym + 27rU/r„, 



where m is unit vector along the line of sight. Here the quantity J 
is proportional to Qhi, hence the visibility correlation scales as the 
square of Qhi- 

We next shift our focus on the redshift space HI power spec- 
trum P H i(k, z) (Bharadwaj & Ali 2004). We calculate P H i(k, z) 
from N-Body simulations and put this in Eqn. 10 to calculate the 
visibility correlation. As simulations have a finite box size, we 
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Table 3. Characteristics of the GMRT, MWA and MWA 5000 that have been used for calculation of the visibility correlation and noise in the visibility 
correlation as well as noise in image are listed here. The values of the system parameters for GMRT and MWA were taken from http://gmrt.ncra.tifr.res.in/ 
and Bowman (2007) respectively. 



Instrument 


Frequency (MHz) 


Ai/ (MHz) B (MHz) 


AU 


00 (degree ) 


T Sya (K) 


A eff (m 2 ) 


GMRT 


610 


1 16 


32 


0.56 


102 


1590 




330 


1 16 


17 


1.03 


106 


1590 




235 


1 16 


12 


1.45 


237 


1590 


MWA 


< 300 


1 32 


2 


9 


68 


5.6 


MWA5000 


<235 


1 32 


1.6 


11.49 


92 


15 




100 



200 



500 



1 000 2000 



U 



Figure 8. This shows the normalized baseline distribution function p(U, v) 
as a function of baseline U for three instruments (GMRT, MWA and 
MWA5000) for frequency v = 330 MHz. 



patch it with a linearly extrapolated power spectrum with a con- 
stant bias at small k. This is especially appropriate for high red- 
shifts where the field of view of the GMRT is larger than the sim- 
ulation box. MWA, of course, has a much larger field of view and 
the statement applies to that as well. 

At sufficiently large scales one can write Pm(k, z) as (Kaiser 
1987) 



P HI (k, z) = b 2 



P(k,z). 



(11) 



Here b is the scale independent bias at large scales and fi{k) = f/b 
where / = is the redshift distortion parameter. P(k, z) is 
matter power spectrum at redshift z and is dominated by the spatial 
distribution of dark matter. 

4.2 System Noise in the Visibility Correlation 

The noise rms. in real part in each visibility for single polarization 
is 



Vrms — 



\/2kBT sys 



(12) 



where T 8ys is the total system temperature, fcs is the Boltzmann 
constant, A e f / is the effective collecting area of each antenna, Av 
is the channel width and At is correlator integration time. This can 
be written in terms of the antenna sensitivity K — A e ff/2 ks as 



V„ 



K^/2AvAt 



(13) 



The noise variance in the visibility correlation is written as (Ali et 
al. 2008) 



2 

°vv = 



' Vrn 



(14) 



where N p is the number of visibility pairs in a particular U bin and 
frequency band B. The HI signal is correlated but the system noise 
is taken to be uncorrelated in visibilities. 

We next consider how N p is calculated. We consider an ele- 
mentary grid with cell of area AU 2 in the u — v plane where the 
signal remains correlated. Total observation time we consider is T. 
We further define p(U, v) (circularly symmetric) to be the baseline 
fraction per unit area per in the baseline range U and U + AU and 
is normalised: J d 2 Up(U, v) — 1. Number of visibility pairs in a 
given U, Av bin in the grid is 



N(N-l) T 
~2 At 



AU 2 p(U, v) 



(15) 



Here N is the total number of antennae Note that the baseline frac- 
tion function p(U, v) is different for different interferometric ar- 
rays and plays an important role in determining the sensitivity of 
an array. This function is also frequency dependent. The signal and 
the system noise is expected to be isotropic and only depends on the 
magnitude of U. We take average over all cells of area AU 2 in the 
circular annulus between U and U + AU. Number of independent 
elementary bins in a circular annulus between U and U + AU is 
2-kUAU/AU 2 . Moreover there are Bj Av independent channels 
for the total bandwidth B. Combining all these we have 



N(N- 



U T -AU 2 p(U,u) 



At 



B 2-kUAU 
A~v AU 2 



(16) 



The noise rms in the visibility correlation can then be written as 

aw = 



2ttN(N - 1) 



K 



1 



tVaUb u - 5 au 15 p{u,v) 

(17) 

The bin size AU is determined by the field of view of the antenna 
and for a Gaussian antenna beam pattern A(6) — e~ e / 8 °, one 
can show that AU = 1/kOq. The equation 17 gives the system 
noise in the visibility correlation for any radio experiment. Note 
that, the system noise rms in the visibility correlation scales as ~ 
l/N(N — 1)T unlike the noise rms in the image which scales as 

~ l/y/N(N-l)T. 

There are two effects that we ignore in the analysis here. 

• At small U, comparable with the inverse of the field of view, 
cosmic variance also limits our ability to measure the visibility cor- 
relations. As the field of view of the GMRT is (much) smaller than 
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the MWA, this concern is more relevant in that case. Cosmic vari- 
ance is subdominant to the system noise as long as we work with 
scales that are much smaller than the field of view, or work with 
U that are much larger than the one that corresponds to the field of 
view. 

• For MWA, the field of view is very large and one needs to 
take the curvature of sky into account (McEwen & Scaife 2008; Ng 
2001, 2005). We ignore this as we are not estimating the signal at 
very large angles. 

We consider three instruments: the GMRT, MWA and a hy- 
pothetical instrument MWA5000. The GMRT has a hybrid distri- 
bution with 14 antennae are distributed randomly in the central 
core of area 1 km x 1 km and 16 antennae are placed along three 
arms (Y shaped). The antennae are parabolic dishes of diameter 
45 m. The MWA will have 500 antennae distributed within a cir- 
cular region of radius 750 m. The antennae are expected to follow 
p an t(r) oc 1 /r 2 distribution with a 20 m flat core. The hypothetical 
instrument MWA5000 is almost similar to the MWA but with 5000 
antennae. We assume antennae are distributed as p an t{r) oc 1/r 2 
within 1000 m radius and with a flat core of radius 80 m (McQuinn 
et al. 2006). The normalized baseline distribution p(U, v) for fre- 
quency v = 330 MHz is presented in the Figure 8 for the three 
instruments (see Datta et al. (2007) for details). Table 3 lists the 
wavebands and other instrument parameters for which analysis is 
done in this paper. We assume the effective antenna collecting area 
(A e f f) to be same as the antenna physical area. We also list param- 
eters used for calculation of noise in visibility correlations. 



4.3 Noise in Images 

Observed visibilities are used to construct an image of the sky for 
each frequency channel. The size of the image is of the order of the 
primary beam 0q and the resolution depends on the largest baselines 
used. Both the size of the image and resolution, or the pixel size 
depend on the frequency of observation in the same manner. 

Unlike visibilities, that are uncorrelated outside of the small 
bins, pixels in an image are not completely uncorrelated. The cor- 
relation of pixels in a raw image is similar to the dirty beam. There- 
fore the analysis of noise in images as a function of scale is a fairly 
complex problem. Scaling of noise with bandwidth is easy as the 
noise for different frequency channels is uncorrelated. raw noise for 
a given pixel can be computed using (assuming two polarizations) 
(Thompson, Moran, & Swenson 2001) 

T 1 

_ -t sys 1 

The symbols in the above equation have the same meaning as de- 
scribed in §4.2. 

We find that for the z — 3.3 window of the GMRT o image ~ 
3.6 /itJy for an integration time of 2000 hours and a bandwidth of 
1 MHz. This has been computed for the central square of the GMRT 
with N = 14. The corresponding number of z = 1.3 is a im age — 
7.8 /iJy for an integration time of 400 hours and a bandwidth of 
1 MHz. 

As a first approximation, we assume that noise in images is 
uncorrelated and that it scales as the square root of the number of 
pixels over which signal is smoothed. This allows us to estimate 
signal to noise ratio for extended structures, but it must be noted 
that this analysis is approximate. 



5 PROSPECTS FOR DETECTION: SIMULATED MAPS 
AND SIGNAL 

In this section we describe results from simulated HI maps as de- 
scribed in the previous sections. 

5.1 Visibility Correlation 

We use the power spectrum of the HI distribution derived from 
our N-Body simulations to compute the visibility correlations. For 
large scales, larger than the simulation size, we patch this with the 
linearly extrapolated power spectrum with a constant bias. This is 
the only input required in Eqn. (10) for calculation of visibility cor- 
relation. 

Visibility correlation for various redshifts is shown in Fig- 
ure 9. We have also plotted the expected system noise in each of 
the panels. We use the the functional form of p(U, v) from Datta 
et al. (2007) (also see Figure 8). The expected visibility correlation 
is shown as a function of U. The expected signal is shown with a 
solid curve and the expected noise in the visibility correlation is 
shows by a dashed curve. This has been shown for three instru- 
ments: GMRT, MWA and a hypothetical instrument MWA5000. 
The noise in each panel has been computed assuming 10 3 hours 
of integration time and visibilities correlated over 1 MHz. We find 
that MWA5000 should detect visibility correlations for baselines 
U < 500 at z ~ 5 and at U < 400 at z ~ 3.7. Higher observation 
time is needed for the MWA to detect the visibility correlations. The 
detection at smaller U should be possible at high significance level. 
Whereas the prospects for detection with the GMRT are encourag- 
ing for 2 = 1.3 where it should be possible to make a detection for 
baselines U < 600. 

It is possible to enhance the signal to noise ratio by combin- 
ing data from nearby bins in U, thus detection may be possible 
in a shorter time scale as well. The GMRT has been in operation 
for more than a decade, and its characteristics are well understood. 
Therefore it is important to make an effort to observe HI fluctua- 
tions at z ~ 1.3 with the instrument. The observations need not 
pertain to the same field as the region observed by the GMRT in a 
single observation is fairly large and we do not expect significant 
fluctuations in the HI power spectrum and hence the visibility cor- 
relation from one field to another. Thus good quality archival data 
for a few fields can be combined, in principle, to detect visibility 
correlations. 

Detection of visibility correlations in turn gives us a measure- 
ment of the power spectrum for the HI distribution (Bharadwaj & 
Sethi 2001). This can be used to constrain various cosmological pa- 
rameters (Bharadwaj, Sethi, & Saini 2009; Visbal, Loeb, & Wyithe 
2008; Loeb & Wyithe 2008). Recently, it has also been pointed out 
that the power spectrum can also be used to estimate the character- 
istic mass of damped lyman alpha systems (DLAS) (Wyithe 2008). 

5.2 Rare Peaks 

We now discuss the possibility of directly detecting rare peaks in 
the HI distribution. This is an interesting possibility for instruments 
like the GMRT and ASKAP that can resolve small angular scales. 
At such scales, non-linear gravitational clustering and the large bias 
for the HI distribution combine to enhance the amplitude of fluc- 
tuations significantly (§3). It is also of interest to see whether this 
enhancement can make detection of rare peaks as easy as statistical 
detection via visibility correlations. 

Figure 10 shows the expected signal to noise ratio (SNR) for 
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Figure 9. This figure shows the expected visibility correlation as a function of U . The expected signal is shown with a solid curve and the expected noise in 
the visibility correlation is shows by a dashed curve. This has been shown for three instruments: GMRT, MWA and a hypothetical instrument MWA5000. The 
noise in each panel has been computed assuming 10 3 hours of integration time and that visibilities are correlated over 1 MHz. We note that the MWA5000 
should detect visibility correlations at U < 500 at z ~ 5 and at U < 400 at z ~ 3.7. The detection at smaller U should be possible at a high significance 
level. Prospects of detection with the GMRT are encouraging for z = 1.3 where it should be possible to make a detection for U < 600. It is possible to 
enhance the signal to noise ratio by combining data from nearby bins in U, thus detection may be possible in a shorter time scale as well. 
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Figure 10. This figure shows the expected signal to noise ratio (SNR) for direct detection of rare peaks in the HI distribution at high redshifts. All the plots 
pertain to the GMRT and we use the expected noise in the image for the central square. The left column is for z ~ 3.3 and the right column is for z ~ 1.3. 
The top row is for variation of the SNR with bandwidth for one pixel, whereas the bottom row shows the variation of SNR with the square root of the solid 
angle over which signal has been smoothed for a given bandwidth. These figures show that a 4 — 5<r detection of a rare peak is possible in 400 hours for 
z ~ 1.3 and in 2000 hours for z ~ 3.3. 



the brightest region in the simulated maps if it is observed using 
the central square of the GMRT. The left column is for z ~ 3.3 and 
the right column is for z ~ 1.3. The top row shows variation of the 
SNR with bandwidth for one pixel, whereas the bottom panel shows 
the variation of SNR with the square root of the solid angle over 
which signal has been smoothed for a given bandwidth. In plots 
of variation of SNR with the bandwidth, the peak is very close to 
the FWHM of spectral lines for large structures (see spectra of rare 
peaks in §7). Although SNR falls off towards smaller bandwidths, 
it does not decline significantly towards larger bandwidths and this 
can be attributed to the presence of nearby lines and is a reflection 
of clustering at small scales. 

These figures show that a 4 — 5a detection of a rare peak is 
possible in 400 hours for z ~ 1.3 and in 2000 hours for z ~ 3.3. 
The time for detection here is comparable to, indeed slightly less 
than the time required for a statistical detection. This is a very ex- 
citing possibility as it allows us to measure the HI fraction of a 
fairly large region, potentially leading to a measurement of Qhi- 
However, unlike statistical detection of the HI distribution, the in- 



tegration time cannot be divided across different fields for a direct 
detection. 

The volume sampled in observation of one field with the 
GMRT at z ~ 3.3 is much larger than the volume of the simu- 
lation used here. Thus there is a strong possibility of finding even 
brighter region in a random pointing at that redshift. On the other 
hand, the volume sampled by the GMRT at z ~ 1.3 is smaller by a 
factor of more than 5 as compared to the volume of the simulation 
used for making mock maps. We have studied the distribution func- 
tion for signal in pixels and find that there are a sufficiently large 
number of bright pixels in the simulated map and a random point- 
ing should have a rare peak comparable to the one discussed here 
within the field of view. This is shown in Figure 1 1 where we have 
shown the number density of pixels in the simulated radio map with 
signal above a given threshold. Here, each pixel corresponds to the 
angular resolution of the central square of the GMRT and a band- 
width of 1 MHz. The volume covered by one GMRT field with a 
bandwidth of 32 MHz should contain more than ten pixels brighter 
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Figure 11. This figure shows the number density of pixels in the simulated 
radio map with signal above a given threshold. Here, each pixel corresponds 
to the angular resolution of the central square of the GMRT and a bandwidth 
of 1 MHz. 

than 10 /xJy. For comparison, the brightest pixel in the simulated 
map gives out twice this flux. 

We do not discuss the possibility of detection of rare peaks 
with the MWA as it has a poorer sensitivity at small scales. 

Note that in our results we have ignored the foregrounds. The 
foregrounds at relevant frequencies will dominate over the cosmo- 
logical HI 21 -cm signal and needs to be subtracted. Line of sight 

modes fcy < 0.065^| (^^p) Mpc" 1 can not be measured 
because of foreground subtraction, B is the bandwidth over which 
foregrounds are subtracted (McQuinn et al. 2006). This would ef- 
fect the baselines U < 80 for a 8 MHz subtraction bandwidth. 
Rare peaks are much smaller objects (< 10 arcmin) than the scales 
where the foreground subtraction would affect and hence we expect 
that our results on rare peaks should not be affected by foregrounds. 



6 DISCUSSION 

Several attempts have been made in recent years to model the 
HI distribution in the post-reionization universe (Scott & Rees 
1990; Subramanian & Padmanabhan 1993; Kumar, Padmanabhan, 
& Subramanian 1995; Bagla, Nath, & Padmanabhan 1997; Bharad- 
waj, Nath, & Sethi 2001; Bharadwaj & Sethi 2001; Bagla & White 
2003; Bharadwaj & Srikant 2004; Bharadwaj & Ali 2005; Loeb & 
Wyithe 2008; Wyithe, Loeb, & Geil 2008; Pritchard & Loeb 2008; 
Wyithe & Brown 2009), with the recent resurgence in interest be- 
ing due to upcoming radio telescopes and arrays. In this paper we 
have revisited the issue using dark matter simulations and an ansatz 
for assigning HI to dark matter haloes. This is amongst the few at- 
tempts made using simulations where the smallest haloes that may 
contain significant amount of HI are resolved, and, the simulations 
are large enough to limit errors due to a finite size of the simulation 
box. 

The prospects for detection of the HI distribution at high red- 



shifts are very promising. The GMRT can be used to measure vis- 
ibility correlations at z ~ 1.3. We can expect the hypothetical 
MWA5000 to measure visibility correlations at z > 3.7 and pro- 
vide an estimate of the power spectrum. The integration time re- 
quired is not very large. This is very encouraging as no other exist- 
ing or upcoming instrument can probe this redshift, though ASKAP 
and MeerKAT will be able to probe the HI distribution at z < 1. 
Measurement of visibility correlations, and hence the HI power 
spectrum can be used to constrain several cosmological parame- 
ters. If we are able to make measurements at several redshifts then 
it becomes possible to constrain models of dark energy. Observa- 
tions of visibility correlations can also be used to put constraints on 
the characteristic mass of damped Lyman alpha systems (Wyithe 
2008). The quadratic dependence of the visibility correlations on 
Qhi can be used to determine the density parameter for HI at high 
redshifts. 

An interesting potential application of the method used here 
is to make predictions for z < 1 where statistical and individual 
detection with the GMRT and ASKAP is likely to be much eas- 
ier. It is obvious that detection of rare peaks of HI with the GMRT 
may be even easier at lower redshifts (z < 0.4). The main issue 
to be addressed at these redshifts is that the observed volume be- 
comes small and hence we can expect significant scatter in the ob- 
served visibility correlations from one field to another, and across 
frequency channels. This is true even for rare peaks and one may 
need to observe several fields before finding a very bright object. Or 
one may need to observe for longer periods of time in order to ob- 
serve a not so rare peak. This is particularly relevant for the GMRT 
where the primary beam covers a solid angle that is nearly 14 times 
smaller than that covered by ASKAP. We propose to study this is- 
sue in a later publication, where we will also discuss prospects of 
observing HI using ASKAP and MeerKAT. We are also working 
on using wavelet based methods for detection of rare peaks in the 
HI distribution. 

In this work we have used a fairly simple HI assignment 
scheme while ignoring the evolution of gas content of haloes. We 
are working on rectifying this shortcoming and expect to use a 
semi-analytical model to improve this aspect of modeling. 

Key conclusions of this paper may be summarized as follows: 

• We find that non-linear gravitational clustering enhances the 
amplitude of perturbations by a significant amount at small scales. 

• HI distribution is strongly biased at high redshifts and this 
enhances the HI power spectrum significantly as compared to the 
dark matter power spectrum. Bias b(k) is scale independent at large 
scales k <JC k n i, where k n i is the scale of non-linearity. 

• Bias decreases sharply from high redshifts towards low red- 
shifts. This leads to a gradual decrease in the brightness tempera- 
ture power spectrum. The change in the amplitude of the brightness 
temperature power spectrum is slow, being less than a factor two 
between z = 5.1 and z = 1.3. 

• Brightest regions in the simulated radio maps are found to be 
extended, with an angular extent much larger than the resolution of 
the GMRT. This can be used to enhance prospects of detection. 

• Spectra in the simulated maps appear to have an FWHM of 
around 1 MHz for z = 1.3 and about half of this for z = 3.34. A 
comparison with earlier work for z = 3.34 shows that the FWHM 
is found to be smaller in the present work. We attribute this to the 
low HI fraction assigned to the most massive haloes, whereas such 
haloes were not excluded in earlier work. 

• The MWA5000 should detect visibility correlations at U < 
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300 at z ~ 5 and at U < 700 at z ~ 3.7. The detection at smaller 
(7 should be possible at a high significance level. 

• Prospects of detection with the GMRT are encouraging for 
z — 1.3 where it should be possible to make a detection for U < 
600. 

• Signal for the brightest source of redshifted 21 cm radiation 
from z = 1.3 appears to be significant, and within reach of an 
instrument like the GMRT. 

• Detection of rare peaks in the HI distribution offers exciting 
possibilities since it does not require a much longer integration time 
as compared to statistical detection. The size of the region repre- 
sented in these rare peaks is fairly large, and it should possible to 
establish the mass contained in these regions using observations in 
other wavebands. The direct detection at z ~ 1.3 requires only a 
few hundred hours of observing time with an existing instrument. 
This may allow us to estimate Qhi in emission, and hence provide 
and independent measurement of the amount of cold gas. 
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